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ABSTRACT 

Context. The M dwarf Gliese 581 has recently been found to harbour two super Earths in addition to an already known close-in Neptune-mass 
planet. Interestingly, these two planets are considered as potentially habitable, and recent theoretical works give further credit to this hypothesis, 
in particular for the outermost planet (Gl 581 d). 

Aims. In this paper, we address the issue of the dynamical stability evolution of this planetary system. This is important because the basic 
stability ensures that a 3-planet model is a physically adequate description of the radial-velocity (RV) data. It is also crucial when considering 
the planets' habitability because the secular evolution of the orbits may regulate their climate, even in the case where the system is stable. 
Methods. We have numerically integrated the planetary system over 10 s yrs, starting from the present fitted solution. We also performed 
additional simulations where i) we vary the inclination of the system relative to the line of sight, ii) assume eccentricities at the upper limit 
of the error bars in the radial velocity fit and where iii) we consider additional (yet undetected) outer planets. We also compute Lyapunov 
exponents to quantify the level of dynamical chaos in the system. 

Results. In all cases, the system appears dynamically stable, even in close to pole-on configurations. The system is actually chaotic, but 
nevertheless stable. The semi-major axes of the planets are extremely stable, and their eccentricities undergo small amplitude variations. The 
addition of potential outer planets does not affect this result. 

Conclusions. Consequently, from the dynamical point-of-view, a 3-planet model is an adequate description of the present RV-data set. Only a 
limited range of inclinations can be excluded for coplanar orbits (i < 10°). The climate on the planets is expected to be secularly stable, thus 
not precluding the development of life. Gl 581 remains the best candidate for a planetary system with planets that potentially bear primitive 
forms of life. 

Key words. Planetary systems - Methods: N-body simulations - Celestial mechanics - Stars: Gliese 581 - Astrobiology - Stars: low-mass, 
brown dwarfs 



1. Introduction 

The M dwarf Gliese 581 has been the subject of a recent in- 
vestigation with the identification of its 3-planet system. One 
of the planets (G1581 b), a Neptune-mass object orbiting the 
star on a 5.4-day or bit, has al ready been known for two years 
dBonfils et alj2005l) . Recentlv. ludrv et all d2007l) have reported 
the discovery of two additional super Earths (G1581 c and d), 
revolving around the star in 12.9 and 83 days (see de t ails in 
Table [B- Considering G1581's luminosity. llJdrv et all d2007h 
inferred the equilibrium temperature of both planets and con- 
clude they may lie within the habitable z one of the star- 
Detailed furth er modeling by Ivon Bloh et alj ( 12007b and 
Selsis et alj d2007l) addre s ses th e habitability of the planets. 
Likewise, Ivon Bloh et al.1 d2007[) find that the greenhouse ef- 
fect increases Gl 58 1 c temperature so much that they no longer 
consider the planet to be in the habitable zone. For similar rea- 
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sons. ISelsis et alj d2007l) find that G1581 c's surface tempera- 
ture is very likely higher than the equilibrium temperature cal- 
culated by lUdrv et al.l (12007b . However, they do not rule out 
habitability for this planet, as a large cloud coverage (> 75%) 
would cool down the planet enough. Conversely, both studies 
agree that the outermost planet (G1581 d) is a good candidate 
for habitability (although close to the outer edge of the hab- 
itable zone) and actually consider it as the better of the two 
candidates. 

An important and unsettled issue about this system con- 
cerns its dynamical behaviour. It is first important to know 
whether the planetary system is dynamically stable and for 
which range of orbital inclinations. If verified, the basic sta- 
bility of the system ensures that the model used (3 planets) is a 
physically adequate description of the observations (the radial- 
velocity measurements). If not verified, the planet detections 
are not necessarily invalidated (as RV periodogrammes clearly 
show that three coherent signals sum up at specific periods). It 
would instead mean that either not enough data were collected 
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Table 1. Orbital parameters of the G1581 planetary system, as derived from the fit of lUdrv et alj d2007l) 



Planet 


Period (days) 


Semi-major axis (AU) 


Eccentricity 


Of (deg) 




f p (JD - 2400000) 


Mass (M ) 


G1581 b 
G1581 c 
G1581 d 


5.36843 ±0.00031 
12.92648 ± 0.00723 
83.22730 ± 0.65845 


0.04061 ±0.16 x 10- 6 
0.07295 ± 2.7 x 10 4 
0.2525 ±0.013 


0.01374 ±0.01405 
0.15926 ±0.05981 
0.12118 ±0.12034 


273.21195 
257.41189 
317.01021 


± 60.54198 
± 24.37209 
±41.79575 


52998.74631 ±0.90393 
52993.40770 ± 1.00130 
52946.80339 ± 11.49888 


15.82 ±0.25 
5.073 ±0.31 
7.804 ± 0.69 



to converge toward the "true" parameters of the system or that 
the model (3 planets) is not complex enough to describe the 
data. Further varying the orbital inclination, we expect to find 
that below a given value - or, equivalently, above given masses 
for the planets - the system becomes unstable. Not valid phys- 
ically, this range of inclinations should be rejected among the 
possible solutions. This partially constrains the sin i degeneracy 
inherent to radial-velocity detections. 

Beyond the basic stability, the secular evolution of the or- 
bits may play an important role regarding planets' habitabil- 
ity. A ll climate calculations ( von Bloh et alj|2007t ISelsis et al.1 
20071) have been done with the currently determined orbits. The 
secular evolution of the orbits has the potential of affecting the 
climate on the planets. A given planet may lie within the hab- 
itable zone but, if subject to significant eccentricity changes, 
it can undergo strong climate variations that eventually pre- 
clude life development. The presently determined eccentrici- 
ties (Table [TJ are small enough to ensure climate stability. But 
one needs to know which maximum values they reach due to 
secular perturbations. 

In the present paper, we numerically investigate the secu- 
lar evolution of the G1581 system, starting from the solution 
of Table Q] In Sect. 2, we study this solution (that we subse- 
quently refer to as the nominal case). In Sect. 3, we perform 
other integrations, assuming different inclinations from edge- 
on and letting the initial eccentricities of the planets reach their 
maximum values within their error bars. In order to quantify 
the dynamical chaos in this system, we compute Lyapunov ex- 
ponents for all these solutions in Sect. 4. In Sect. 5, we investi- 
gate the perturbing action of potentially additional outer plan- 
ets that have not been detected yet, provided their contribution 
to the radial velocity signal is small enough compared to the 
residuals of the 3-planets fit. Our conclusions are presented in 
Sect. 6. 

2. The nominal case 

The best 3-planet orbital fit for G1581 is explained in Table Q] 
This solution with the assumptions of coplanarity and sin i = 1 
(i = 90°, an edge-one system) will constitute our nominal 
case. We numerically integrate this system taking 0.31 M Q 
for the mass of G1581. The integratio n is performed usin g 



the symplectic N-body code SyMBA (IDuncan et al.l [1998), 



which handles close encounters. The initial timestep is fixed 
to 2 x 10~ 4 yr = 0.18 day, i.e. 1/30 of the smallest orbital 
period. Symplectic integration schemes usually ensure en- 
ergy conservation with 10~ 6 relative accuracy as soon as the 
timestep is taken to ~ 1/20 of the smallest orbital period 
dLevison & Duncanlll994l) . The integration is carried out over 
10 8 yr. On more limited timespans, we checked that taking a 
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Fig. 1. Fractional errors on the total energy E with respect to the 
initial one Eq (black curve), and on the total angular momentum 
H with respect to the initial one Ho (grey curve), as a function 
of time over the 10 8 yr integration 

Table 2. Precession frequencies for the nominal solution, as 
computed from the linear secular theory 



Name 


Frequency ("/yr) 


Period (yr) 


gi 


3300.9 


392.62 


82 


539.04 


2404.2 


83 


38.199 


33945. 



significantly smaller timestep does not change the result. In 
Fig.[T] we display the fractional errors on the total energy and 
angular momentum over the 10 8 yr integration. The energy is 
preserved to less than 10~ 7 relative accuracy. Hence we are con- 
fident in our integration. Figure [2] shows the first 10 4 years of 
the integration. We see that the secular variations of the 3 plan- 
etary orbits are very regular. The eccentricities undergo quasi- 
periodic modulations, while the longitudes of periastra precess 
regularly. This solution is in fact very close to the one we can 
compute with a linear sec ular theory (Laplace - Lagrange), 
such as described the one by lBretagnorJ(ll974ll990l) . In the lin- 
ear approximation, the secular evolution of the eccentricity vec- 
tors of the 3 planets is a combination of sine and cosine terms 
with 3 characteristic frequencies (gi, i — 1,2,3) that are listed 
in Table|2] These frequencies are obtained by solving the linear 
secular equations for their eigenvalues. We obviously see these 
characteristic frequencies in Fig. [2] An interesting outcome is 
that these precession frequencies are much higher than in the 
Solar System, which do not exceed 25"/yr, and g\ is basically 
the main precession frequency of the periastron of G1581 b 
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Fig. 2. The 10 4 first years of the integration of the nominal solution. The upper plots show the temporal evolution of the eccen- 
tricity of the three planets G1581 b,c, and d, from left to right, respectively; the lower plots are the same for the longitude of 
periastron m 



Table 3. Variation ranges for some orbital parameters fo the 3 planets over the 10 yr integration. 



Planet 


Semi-major axis 


Eccentricity 


Periastron 


Apoastron 




(AU) 




(AU) 


(AU) 


G1581 b 


0.040609 - 0.046185 


0.01 - 0.095 


0.0368 - 0.0402 


0.041 -0.0445 


G1581 c 


0.072885 - 0.073 


0.07-0.16 


0.0614 - 0.0678 


0.078 - 00846 


G1581 d 


0.2522 - 0.2528 


0.1200 - 0.1246 


0.2207 - 0.2227 


0.2823 - 0.2843 



and G1581 c. This is obviously due to the much smaller size 
of the system. Given the error bar on the fits of the arguments 
of periastra {a>\, 0)2), the secular motion should be detectable 
within ~ 30 years, and probably less if the orbital fits get more 
constrained in the near future thanks to further monitoring. 

In Table|3] we list the maximum evolution ranges for the or- 
bital elements of the three planets. The semi-major axes are ex- 
tremely stable, revealing a regular dynamics out of any mean- 
motion resonance configuration. The evolution ranges of the 
eccentricities are narrow, so that we may claim than the system 
is stable with a high level of confidence. While the time span 
of the integration is 10 8 yr, most of the characteristic features 
of the secular evolution of the orbital parameters occur on a 
10 4 yr-time scale. Therefore, even if the star is believed to be 
older than 2 x 10 9 yr, the current integration clearly explores all 
the dynamical possible outcomes of the system. Actually, due 
to the short orbital periods of the planets (and to the high pre- 
cession frequencies), integrating the G1581 over 10 s yris basi- 
cally equivalent to integrating the Solar System over ~ 100 Gyr 
I 

Interestingly, the present-day eccentricity of G1581 c 
roughly corresponds to its maximum values along its secular 
evolution, and the eccentricity of G1581 d only has small vari- 
ations. Hence we expect the climate of both outer planets to be 
secularly stable. 



3. Other solutions 

The nominal solution corresponds to an inclination i - 90° (so 
the lowest possible planetary masses) and to the orbital parame- 
ters of the discovery paper (TableQ]). Lower inclinations and/or 
parameter's values slightly outside the best solution may lead 
to different dynamical behaviours that are worth investigating. 

In a first set of additional simulations, we assume various 
inclinations ranging from (pole on) to 90° (edge on), but still 
holding the initial eccentricities to their nominal values. The 
mass of each planet is augmented by a factor 1 / sin i with re- 
spect to the nominal case. In a second set of simulations we 
assume different inclinations and, moreover investigate the im- 
pact of eccentricities larger than in the nominal case (as less 
stability is only expected if the eccentricity is larger). For that 
set, we take the initial eccentricities for the three planets at the 
upper limit of their error bars (we add 1 <x to the eccentricities) 
For both sets of integrations, we plot the width of the evolution 
ranges obtained over the 10 5 yr integration for both the semi- 
major axis and the eccentricities of the three planets (Fig. [3). 

As can be seen from Fig. [3] when the inclination decreases, 
the dynamical interactions increase accordingly and we ex- 
pect the system to become unstable below a given inclination. 
As for the nominal case, the integrations are carried out over 
10 5 yr. They naturally show that both the semi-major axis and 
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Fig. 3. Stability of the three-planet system in various configurations. The maximum variation range for the semi-major axes 
(upper plots) and for the eccentricities (lower plots) is displayed as a function of the assumed viewing inclination of the system 
with respect to pole-on. Each cross corresponds to a single simulation. The left plots correspond to simulations with the nominal 
eccentricities as initial conditions, and the right plots to simulations with eccentricities increased by lcr relative to the error bars 
given in Table [T] 



the eccentricity take a wider range of values than in the nom- 
inal case with decreasing inclinations. In the first set of inte- 
grations (nominal initial eccentricities), the system nonetheless 
remains stable down to i = 10°. Almost pole-on configurations 
(i < 10°) are unstable and should be rejected from possible so- 
lutions. Although, such low inclinations are very improbable, 
and from the statistical point-of-view, the actual masses of the 
planets are probably close to those listed in TableQ] 



In the second set of simulations (lcr augmented initial ec- 
centricities), the dynamical interactions are slightly enhanced 
and the semi-major axis and the eccentricity take a wider range 
of values than for the first set of additional simulations. The 
system is therefore found unstable below larger inclinations 
(< 20°). 



In all cases, the instability appears very unlikely. If we 
assume that the rotation axis of the system is randomly dis- 
tributed in space, i > 20° occurs with a probability of 0.94. 
In conclusion, irrespective of its actual inclination, the G1581 
planetary system is very probably stable. 
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Fig. 4. Progress of the calculation of the partial Lyapunov ex- 
ponents as a function of the integration time, in the nominal 
case for G1581 (zero inclination, nominal eccentricities). At 
t — 10 4 yr, the three exponents have stabilised. 
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Fig. 5. Lyapunov exponents as a function of the inclination i, computed for all the simulations described in Fig. [3] Left plot 
simulations with nominal eccentricities; Right plot : simulations with lcr increased eccentricities. 



4. Lyapunov exponents 

Looking for variation ranges for the orbital elements is a ba- 
sic tool for investigating the stability of a planetary system. 
A more sophisticated way for quantifying chaos is to com- 
pute Lyapunov exponents. For all simulations decribed above, 
we compute the maximum Lyapunov ex ponents (MLE) for the 
three pla nets, following the technique bv lBenettin et al.1 dl978l) 
(see also iMorbidellil 120021) . The exponents are computed by 
integrating fictitious bodies having initial conditions that are 
very close to those of the planets and estimating their exponen- 
tial diverge rate. We coupled this algorithm with the SyMBA 
integrator. 

When we start integrating the bodies with initial coordinate 
vector po (holding for the positions and velocities of all the 
bodies), we also integrate another system of bodies with iden- 
tical masses, but with initial coordinate vector po + 6po, where 
II^Poll ^ llPoH- After a fixed normalization time f n0 rm, we com- 
pute the error vector 6p(t aorm ) as the difference at t — f norm (af- 
ter integration) between the coordinate vector of the fictitious 
bodies and of the regular bodies. We then compute 

ll<5/»( f norm)ll 



Si = 



and 



dpi = 



8p(tnorm) 



(1) 



\\dpo\\ Si 

We then use dpi as a new initial error vector for the ficti- 
tious bodies relative to the coordinates vector of the regular 
bodies at t norm , and we iterate the above process. Each f norm , 
the error vector is renormalized this way, and we obtain a 
sequen ce si, S2, 



of renormalization factors. Benettin et al 



(1978) proved that the MLE X. can be computed as 



£ = 



lim V* Si = lim L„ 



(2) 



i=i 



The result is independent of the choice of f n orm, provided it 
is chosen small enough to avoid too large an exponential di- 
vergence. During the integration, we compute log L„ for every 
£ norm , and we try to derive an asymptotic behaviour. Two cases 
can occur: i) log L„ converges towards a finite limit. Then the 
system is chaotic and we have reached X. > within our in- 
tegration time, ii) Or logL„ keeps decreasing monotonically 



at the end of the integration. This means that the orbit is very 
probably regular or, more precisely, less chaotic than a given 
level that depends on the integration time. 

In practice we compute MLEs for all the individual bod- 
ies in the integration: for each orbit we compute the associated 
(partial) error vector, and perform individual renormalizations. 
Each body has its own sequence of s, 's. Note that more than the 
absolute values of the MLE derived, the comparison between 
the values derived for the individual bodies and between dif- 
ferent integrations is more relevant. This shows clearly which 
orbits are the most chaotic. 

In our example applications, we integrate over 10 4 yr to 
compute the MLEs. f norm has been fixed to 0.02 yr, i.e., 100 
times the timestep. The progress of the computation of the 
MLEs in the nominal case for the three planets is plotted as 
a function of the integration time in Fig. [4] We see that they 
have all converged towards finite limits at t — 10 4 yr, showing 
that the orbits are actually chaotic. 

The global result of MLE calculation is shown in Fig. [5] 
where we have computed the MLEs for the three planets for 
all the simulations described in Fig.|3](stopping at t — 10 4 yr). 
In all cases, we obtain non-zero exponents, showing that the 
system is actually chaotic. 

We see that the MLEs slowly increase with decreasing in- 
clinations, showing as expected that solutions at smaller incli- 
nations are more chaotic, due to higher planetary masses. We 
nevertheless note that the variation is small except for i < 20°. 
The system is not much more chaotic at i = 20° than at i = 90°. 
This confirms that there is no real dynamical constraint on the 
inclination. We also note that solutions with lcr increased ec- 
centricities are not more chaotic than those with nominal ec- 
centricities. As a result the dynamical stability does not put 
any additional constraint on the planet eccentricities other than 
those derived from the radial velocity analysis. 

From Fig. [5] it also becomes clear that the two inner plan- 
ets (G1581 b and c) are much more chaotic than the outer one 
(Gl 58 1 d) (the exponent is smaller). In fact, the two inner plan- 
ets are significantly chaotic. This does not prevent them from 
being stable. Actually, chaos does not necessarily mean insta- 
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Table 4. Semi-major axis and eccentricity variation ranges for 
the 3 known planet of Gl 58 1 , plus additional outer planets (see 
text), computed over 10 5 yr integrations 



Planet 



Semi-major axis (AU) 



Eccentricity 



With a 1 Mi planet at 5 AU : 

G1581b 0.04060120-0.04061199 
G1581c 0.07286312-0.0729903 
G1581d 0.25222823 -0.25281971 



0.00367502 - 0.0895939 
0.09259277 -0.164875 
0.11884516-0.1228305 



With a 29.6 M e planet at 5 AU : 

G1581b 0.04060139-0.04061188 
Gl 58 1 c 0.07286522 - 0.0729892 
G1581 d 0.25222915 -0.2528203 



0.00347993 -0.0895106 
0.09258738 -0.1647922 
0.11887163 -0.1226443 



With a 26.5 M e planet at 4 AU : 
G1581b 0.04060141 -0.04061174 
Gl 58 1 c 0.07286447 - 0.07298892 
Gl 58 1 d 0.25222689 - 0.25282013 



0.00354749-0.08949917 
0.09258812-0.16479787 
0.11884819-0.1226468 



With a 22.9 M e planet at 3 AU : 

G1581 b 0.04060137 -0.04061170 
Gl 58 1 c 0.07286369 - 0.07298806 
Gl 58 1 d 0.25222659 - 0.2528 195 1 



0.00372438 -0.0895144 
0.09258945 -0.16479738 
0.11884226-0.12267108 



With a 18.7 M e planet at 2 AU : 

G1581b 0.04060140-0.04061177 
G1581c 0.07286632-0.07298852 
G1581d 0.25222099-0.25281587 



0.00405747 - 0.08955927 
0.09259499-0.16492574 
0.11886586-0.12266511 



With a 13.2 M e planet at 1 AU : 

G1581b 0.04060131 -0.04061176 
Gl 58 1 c 0.07286289 - 0.07298820 
Gl 58 1 d 0.25222239 - 0.2528001 1 



0.00475611 -0.08969945 
0.09258994-0.16491917 
0.11882570-0.12246267 



bility. The Solar System is known to be chaotic (but on longer 
timescales), yet it is nevertheless stable. 



5. Other planets 

Our simulations were made with the three known planets orbit- 
ing G1581. However, the system may harbour additional, un- 
known planets. The presence of these planets may affect the 
stability of the whole system. There are upper limits to the pres- 
ence of additional (mainly outer) planets. The maximum am- 
plitude of the residuals in the 3-planet fits of lUdrv et al.l Q2007) 
is +2.1 ms~'. Any additional planet should not generate a ra- 
dial velocity with a larger amplitude, otherwise it would have 
already been detected. Assuming i = 90° and a circular orbit, 
this puts severe constraints on the mass m and distance d of the 
unseen planet. We derive 



1M„ 



< 13.227 x 



1 AU 



(3) 



This constraint holds if the unseen planet generates full- 
amplitude variations withi n the timespan o f the radial veloc- 
ity data, i.e., ~ 1000 days dUdrv et alj|2007l) . This means that 



the orbital period of the unseen planet must not exceed -twice 



this time span to account for this constraint, i.e. an orbital dis- 
tance d < 5.5 AU. For more distant planets, the constraint is 
much weaker. In fact, this upper limit is probably already too 
large. Dynamically speaking, we do no expect any hypothetical 
planet orbiting at 5.5 AU to significantly affect the dynamics of 
the inner system located inside 0.25 AU. We are thus confident 
in the conclusions we derive below, as all potentially destabi- 
lizing configurations have been explored. 

We thus performed new simulations, each of them with the 
nominal conditions, but to which we add an additional planet 
orbiting the star on a circular orbit at an arbitrary distance d, 
and with the maximum mass allowed by Eq.©. All the in- 
tegrations were carried out over 10 5 yr. We did 5 simulations 
with d — 5, 4, 3, 2, and 1 AU. This gives masses of 29.6, 26.5, 
22.9, 18.7, and 13.2 M e , respectively. We also added a simula- 
tion with a 1 Jupiter mass (Mj) planet orbiting the star at 5 AU, 
as the constraint (01 is less severe at this distance. Note that 
this case is by far the worst possible disturbing configuration 
that is still compatible with the constraints. More distant com- 
panions, even massive, are less destabilizing. In a first-order 
approximation, the perturbing effect of a distant planet of mass 
m orbiting at distance d on an inner planet orbiting at distance 
r scales as the tidal stripping effect on the orbit, i.e. oc mr/d 3 . 
Hence a 1 Mj planet at 5 AU is as disturbing as a 8 Mj planet 
at 10AU and a 32 Mj brown dwarf at 20 AU. The AO surveys 
would likely have already detected such a massive companion. 

In all cases, the whole system appears just as stable as with- 
out any additional planet. The result concerning the stability 
is summarised in Table [4] where we give the semi-major axis 
and eccentricity variation ranges for the three known planets. 
The results are very similar among the 6 different integrations, 
even in the case of a Jovian planet, showing that the additional 
planet has little influence on the stability of the inner system. 
Moreover, Table |4] is easily compared to Table [3] The varation 
ranges are very similar. Therefore, we may stress that any ad- 
ditional outer planet that fits into the constraint of the radial 
velocity residuals does not affect the stability of the 3-planet 
system. Note that the maximum eccentricity values in Table [4] 
are actually slightly lower than those in Table[3] This could ap- 
pear surprising, since the integration in Table [3] is made with- 
out any additional perturber. Recall, however, that it extends 
over 10 s yr instead of 10 5 for those in Table|4] This shows con- 
versely that if we were entending the integrations of Table [4] 
up to 10 8 yr, we should expect slightly wider variation ranges. 
The basic conclusion nevertheless remains: the stability of the 
system is not affected. 

6. Discussion 

We have computed the secular evolution of the G1581 plane- 
tary system in various possible configurations. The main con- 
clusion is that the system is almost always stable. It is stable for 
inclinations as low as ~ 20° and even if the initial eccentricities 
are augmented by their 1-cr error bars. 

As expected for any planetary system with regular dynam- 
ics, the semi-major axes vary very little and the three planets 
are expected to remain at their current location with respect to 
the star. Meanwhile, the eccentricities of the two outer plan- 
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ets (both considered for habitability) reach values that are sig- 
nificantly above the Earth's value. Concerning G1581 c, the 
present-day eccentricity is close to its maximum value. This 
planet is not expected to get much farther away from its parent 
star and, to maintain a surface temperature cool enough to al- 
low the presence of liquid water, a high water-cloud coverage 
(~75%) would be required at any time. Regarding Gl 58 1 d, the 
nominal eccentricity is non negligible (~ 0.12) and also found 
to be very stable. It is significantly above the maximum value 
reached by the Earth throughout its secular evolution (~ 0.06, 
see e.g. lLaskanll988b and corresponds to a 24% variation of 
the radiation flux received from the star between apoastron and 
periastron. The anomalistic season effects should therefore be 
strong, if not damped by the short orbital period (83 days). If 
we compare the periastron and apoas tron values of G15 81 d 
to the habitabl e zone calculations by ISelsis et al.l ([2007) and 
von Bloh et alJ (120071) . we see that Gl 58 1 d is outside the habit- 
able z one at apoast r on bu t well inside at periastron. As pointed 
out by lSelsis et al.l (120071) . the average stellar flux received by 
an eccentric orbit is enhanced by a factor 1 / Vl - e 2 with re- 
spect to a circular orbit with the same semi-major axis. This 
can help maintaining G1581 d in the habitable zone. What we 
show here is that this effect is secularly permanent. 

Now, if the obliquity of the rotation axis of this planet is 
non-zero, this should combine with the obliquity's seasonal ef- 
fect and lead to climate differences between the hemispheres 
of this planet, much like Mars pres ently. The o bliqui ty of 
G1581 d is of course unknown, but ISelsis et aD (12007) and 
von Bloh et al.l d2007l) agree in claiming that, given the esti- 
mated age of the star (>2 Gyrs), the rotation of Gl 58 1 d should 
already be tidally locked with the orbital motion. In that case, 
we would expect the obliquity to have been set to zero by tidal 
effects, and there should instead be climate differences between 
the night and day hemispheres. This could help in maintain- 
ing the day hemisphere habitable. Sels is et al. ( 20071) also show 
that tidal locking does not contradict the non-zero eccentricity 
of the orbit. Tides usually tend to both synchronize the rota- 
tion and circularize the orbit. The circularization time is almost 
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always l onger than the sync hronization time (lHutl ll981). For 
Gl 581 d. lSelsis et al.l d2007l) estimate the synchronization time 
to lOMyrs and the circularization time to 10 Gyrs, i.e. well 
above the present age of the system. 
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